DTI

Data, EDA

plot, mean

library(refund)
data(DTI); attach(DTI)
names(DTI)

## [1] "ID"         "visit"      "visit.time" "Nscans"     "case"      
## [6] "sex"        "pasat"      "cca"        "rcst"

DTI.complete <- subset(DTI, complete.cases(DTI)) #NA 처리
DTI.baseline <- subset(DTI.complete, visit == 1 & case == 1)
tract <- 1:93
n <- length(unique(DTI.baseline$ID)) #66
dim(DTI.baseline$cca) #66 93

## [1] 66 93

tract2 <- 1:55
dim(DTI.baseline$rcst) #66 55

## [1] 66 55

matplot(tract, t(DTI.baseline$cca), 
        type='l', lty=1, col=rainbow(n),
        main = "Diffusion Tensor Imaging1 : CCA",
        xlab="tract", ylab="Fractional anisotropy (FA)")

matplot(tract2, t(DTI.baseline$rcst), 
        type='l', lty=1, col=rainbow(n),
        main = "Diffusion Tensor Imaging2 : rcst",
        xlab="tract2", ylab="Fractional anisotropy (FA)")

set.seed(245)
n.crv <- 5
sel.crv <- sample(1:n, size=n.crv, replace = FALSE)

matplot(tract, t(DTI.baseline$cca), 
        type='l', lty=1, col="light grey",
        main = "Diffusion Tensor Imaging1 : CCA",
        xlab="tract", ylab="Fractional anisotropy (FA)")

matlines(tract, t(DTI.baseline$cca[sel.crv,]), 
         type='l', lty=1, col=rainbow(n)[sel.crv])

matplot(tract, t(DTI.baseline$cca), type='l', lty=1, col='lightgrey',
        main = "Diffusion Tensor Imaging1 : CCA",
        xlab="tract", ylab="Fractional anisotropy (FA)")
sampleMean <- colMeans(DTI.baseline$cca)
lines(tract, sampleMean, lty=2, lwd=2, col='red')

matplot(tract2, t(DTI.baseline$rcst), type='l', lty=1, col='lightgrey',
        main = "Diffusion Tensor Imaging2 : rcst",
        xlab="tract", ylab="Fractional anisotropy (FA)")
sampleMean <- colMeans(DTI.baseline$rcst)
lines(tract2, sampleMean, lty=2, lwd=2, col='red')

## 2D rainbow plot 
clrs <- rev(colorRampPalette(c("blue", "green", "yellow", "red"))(40))    
colfct <- as.numeric(cut(DTI.baseline$pasat, 40))
matplot(tract, t(DTI.baseline$cca), 
        type = 'l', col = clrs[colfct], lty = 1, 
        ylab = "Fractional Anisotropy",
        xlab = "Distance Along Tract", ylim = range (.2, .85))

## 3D rainbow plot
par(mar=c(1,1,0,0), cex.axis=1, cex.lab=1)
clrs <- rev(colorRampPalette(c("blue", "green", "yellow", "red"))(40))    
?rev
proj = persp(x = tract, y = seq(min(DTI.baseline$pasat), max(DTI.baseline$pasat), l=length(DTI.baseline$pasat)),  z=t(DTI.baseline$cca),
             xlab="tract", ylab="PASAT", zlab="FA", col=NA, border=NA,
             ticktype = "detailed", axes=TRUE, theta=30, phi=30)

o <- rev(order(DTI.baseline$pasat))
for(i in o){
  lines(trans3d(x = tract, y=rep(DTI.baseline$pasat[i], ncol(DTI.baseline$cca)),  z=DTI.baseline$cca[i,], pmat=proj), col=clrs[colfct[i]])
}

cov

library(fields)

## Loading required package: spam

## Loading required package: dotCall64

## Loading required package: grid

## Spam version 2.6-0 (2020-12-14) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction 
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.

## 
## Attaching package: 'spam'

## The following objects are masked from 'package:base':
## 
##     backsolve, forwardsolve

## See https://github.com/NCAR/Fields for
##  an extensive vignette, other supplements and source code

sampleCov <- cov(DTI.baseline$cca)
image.plot(tract, tract, sampleCov, main='sample covariance1 of FA: CCA')

sampleCov1 <- cov(DTI.baseline$rcst)
image.plot(tract2, tract2, sampleCov1, main='sample covariance2 of FA: rcst')

dim(cov(DTI.baseline$cca))

## [1] 93 93

tick <- seq(1,93, by = 10)
persp(tract[tick], tract[tick], sampleCov[tick,tick],
            xlab="track", ylab="track", zlab="", 
      main="Sample Covariance of FA: CCA",
      theta = 30, phi = 30, col='light grey', shade=0.1)

tick2 <- seq(1,55, by = 10)
persp(tract2[tick2], tract2[tick2], sampleCov1[tick2,tick2],
            xlab="track", ylab="track", zlab="", 
      main="Sample Covariance of FA: rcst",
      theta = 30, phi = 30, col='light grey', shade=0.1)

library(lattice)
x <- rep(tract[tick], length(tract[tick]))
y <- rep(tract[tick], each = length(tract[tick]))
z <- as.vector(sampleCov[tick, tick])
dat <- data.frame(x=x, y=y, z=z)
wireframe(z ~ x*y, data =dat, 
  xlab="tract", ylab="tract", zlab="",
  main="Sample Covariance of FA: CCA",
  drape = TRUE,
  colorkey = TRUE,
  screen = list(z = -60, x = -60)
)

sampleCor <- cor(DTI.baseline$cca)
image.plot(tract, tract, sampleCor, main='sample correlation of FA: CCA')

sampleCor2 <- cor(DTI.baseline$rcst)
image.plot(tract2, tract2, sampleCor2, main='sample correlation of FA: rcst')